Paso 5- Trayectorias de hospitalización y mortalidad con foco en condiciones vinculadas a trastornos de salud mental y consumo de sustancias posterior a un primer ingreso por alguno de estos trastornos, en usuarios/as jóvenes y adultos emergentes de población general y pertenecientes a pueblos originarios, 2018-2021, Chile
Representar las mejores opciones de agrupamiento, junto con su relación con otras variables (para poster)
Author
Andrés González Santa Cruz
Published
October 2, 2024
Configurar
Code
# remover objetos y memoria utilizadarm(list=ls());gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 596638 31.9 1347497 72 686445 36.7
Vcells 1122406 8.6 8388608 64 1876197 14.4
#elegir repositorioif(Sys.info()["sysname"]=="Windows"){options(repos =c(CRAN ="https://cran.dcc.uchile.cl/"))}options(install.packages.check.source ="yes") # Chequea la fuente de los paquetes#borrar caché#system("fc-cache -f -v")if(!require(pacman)){install.packages("pacman");require(pacman)}pacman::p_unlock(lib.loc =.libPaths()) #para no tener problemas reinstalando paquetesif(Sys.info()["sysname"]=="Windows"){if (getRversion() !="4.4.0") { stop("Requiere versión de R 4.4.0. Actual: ", getRversion()) }}if(!require(job)){install.packages("job");require(job)}if(!require(kableExtra)){install.packages("kableExtra");require(kableExtra)}if(!require(tidyverse)){install.packages("tidyverse");require(tidyverse)}if(!require(cluster)){install.packages("cluster"); require(cluster)}if(!require(WeightedCluster)){install.packages("WeightedCluster"); require(WeightedCluster)}if(!require(devtools)){install.packages("devtools"); require(devtools)}if(!require(TraMineR)){install.packages("TraMineR"); require(TraMineR)}if(!require(TraMineRextras)){install.packages("TraMineRextras"); require(TraMineRextras)}if(!require(NbClust)){install.packages("NbClust"); require(NbClust)}if(!require(haven)){install.packages("haven"); require(haven)}if(!require(ggseqplot)){install.packages("ggseqplot"); require(ggseqplot)}if(!require(gridExtra)){install.packages("gridExtra"); require(gridExtra)}if(!require(Tmisc)){install.packages("Tmisc"); require(Tmisc)}if(!require(factoextra)){install.packages("factoextra"); require(factoextra)}if(!require(stargazer)){install.packages("stargazer"); require(stargazer)}if(!require(gtsummary)){install.packages("gtsummary"); require(gtsummary)}if(!require(lmtest)){install.packages("lmtest"); require(lmtest)}if(!require(emmeans)){install.packages("emmeans"); require(emmeans)}if(!require(fpp2)){install.packages("fpp2"); require(fpp2)}if(!require(purrr)){install.packages("purrr"); require(purrr)}if(!require(forecast)){install.packages("forecast"); require(forecast)}if(!require(magrittr)){install.packages("magrittr"); require(magrittr)}if(!require(foreach)){install.packages("foreach"); require(foreach)}if(!require(doParallel)){install.packages("doParallel"); require(doParallel)}if(!require(progressr)){install.packages("progressr"); require(progressr)}if(!require(chisq.posthoc.test)){devtools::install_github("ebbertd/chisq.posthoc.test")}if(!require(rstatix)){install.packages("rstatix"); require(rstatix)}seq_mean_t_dos_grupos <-function(bd =NULL, group1, group2) {# Agrupar por ambas variables resultados <-by(bd, list(group1, group2), seqmeant)# Obtener todas las combinaciones posibles de los grupos combinaciones <-expand.grid(group1 =unique(group1), group2 =unique(group2), stringsAsFactors =FALSE)# Extraer los resultados y asociarlos con las combinaciones resultados_df <-do.call(rbind, lapply(seq_along(resultados), function(i) { group_name1 <-attr(resultados, "dimnames")[[1]][i] group_name2 <-attr(resultados, "dimnames")[[2]][i]data.frame(factor_inclusivo_1 = group_name1, factor_inclusivo_2 = group_name2, Mean = resultados[[i]]) }))# Unir los resultados con las combinaciones para rellenar los valores faltantes final_df <-merge(combinaciones, resultados_df, by.x =c("group1", "group2"), by.y =c("factor_inclusivo_1", "factor_inclusivo_2"), all.x =TRUE)return(final_df)}multinom_pivot_wider <-function(x) {# check inputs match expectatations# create tibble of results df <- tibble::tibble(outcome_level =unique(x$table_body$groupname_col)) df$tbl <- purrr::map( df$outcome_level,function(lvl) { gtsummary::modify_table_body( x, ~dplyr::filter(.x, .data$groupname_col %in% lvl) %>% dplyr::ungroup() %>% dplyr::select(-.data$groupname_col) ) } )tbl_merge(df$tbl, tab_spanner =paste0("**", df$outcome_level, "**"))}best_subset_multinom <-function(y, x.vars, data) {# y Nombre de la variable dependiente (cadena de texto)# x.vars Vector de nombres de predictores (caracter)# data Dataframe con los datos de entrenamiento# Cargar las librerías necesariasrequire(dplyr)require(purrr)require(tidyr)require(nnet)require(MASS)# Generar todas las combinaciones posibles de predictores predictors_list <-lapply(1:length(x.vars), function(i) {combn(x.vars, i, simplify =FALSE) }) %>%unlist(recursive =FALSE)# Inicializar una lista para almacenar los resultados results <-list()# Iterar sobre cada combinación de predictoresfor (i inseq_along(predictors_list)) { predictors <- predictors_list[[i]] formula <-as.formula(paste(y, "~", paste(predictors, collapse ="+")))# Ajustar el modelo multinomial model <-tryCatch( nnet::multinom(formula, data = data, trace =FALSE),error =function(e) NULL )# Si el modelo se ajustó correctamente, almacenar los resultadosif (!is.null(model)) {# Extraer el AIC del modelo aic <-AIC(model)# Almacenar la información en una lista results[[length(results) +1]] <-list(predictors = predictors,model = model,AIC = aic ) } }# Convertir la lista de resultados en un dataframe results_df <- results %>% purrr::map_df(function(res) {data.frame(predictors =paste(res$predictors, collapse ="+"),AIC = res$AIC,stringsAsFactors =FALSE ) })# Ordenar los modelos por AIC de menor a mayor results_df <- results_df %>%arrange(AIC)return(results_df)}best_subset_multinom_interactions <-function(y, x.vars, data) {# y Nombre de la variable dependiente (cadena de texto)# x.vars Vector de nombres de predictores (caracter)# data Dataframe con los datos de entrenamiento# Cargar las librerías necesariasrequire(dplyr)require(purrr)require(tidyr)require(nnet)require(MASS)# Generar todas las combinaciones posibles de predictores (efectos principales) main_effects_list <-lapply(1:length(x.vars), function(i) {combn(x.vars, i, simplify =FALSE) }) %>%unlist(recursive =FALSE)# Inicializar una lista para almacenar los resultados results <-list()# Iterar sobre cada combinación de efectos principalesfor (main_effects in main_effects_list) {# Generar términos de interacción de hasta 3 variables interaction_terms <-list()# Para interacciones de 2 variablesif (length(main_effects) >=2) { interaction_terms_2way <-combn(main_effects, 2, function(x) paste(x, collapse =":")) interaction_terms <-c(interaction_terms, interaction_terms_2way) }# Para interacciones de 3 variablesif (length(main_effects) >=3) { interaction_terms_3way <-combn(main_effects, 3, function(x) paste(x, collapse =":")) interaction_terms <-c(interaction_terms, interaction_terms_3way) }# Combinar efectos principales e interacciones all_terms <-c(main_effects, interaction_terms)# Generar todas las combinaciones posibles de términos (incluyendo interacciones)# Solo se incluyen interacciones si sus efectos principales están presentes term_combinations <-list()# Obtener todos los subconjuntos de efectos principales main_effects_subsets <-lapply(1:length(main_effects), function(i) {combn(main_effects, i, simplify =FALSE) }) %>%unlist(recursive =FALSE)# Para cada subconjunto de efectos principalesfor (me in main_effects_subsets) {# Iniciar con los efectos principales terms <- me# Incluir interacciones solo si todos sus efectos principales están incluidos possible_interactions <- interaction_terms[sapply(interaction_terms, function(x) { vars_in_interaction <-unlist(strsplit(x, ":"))all(vars_in_interaction %in% me) }) ]# Generar todas las combinaciones de interacciones para incluir interaction_subsets <-list(NULL)if (length(possible_interactions) >0) { interaction_subsets <-lapply(1:length(possible_interactions), function(i) {combn(possible_interactions, i, simplify =FALSE) }) %>%unlist(recursive =FALSE) }# Para cada combinación de interacciones, crear el conjunto completo de términosfor (ints in interaction_subsets) {if (is.null(ints)) { full_terms <- terms } else { full_terms <-c(terms, ints) }# Añadir a la lista de combinaciones de términos term_combinations <-append(term_combinations, list(full_terms)) } }# Ajustar modelos para cada combinación de términosfor (terms in term_combinations) { formula <-as.formula(paste(y, "~", paste(terms, collapse ="+")))# Ajustar el modelo multinomial model <-tryCatch( nnet::multinom(formula, data = data, trace =FALSE),error =function(e) NULL,warning =function(w) NULL )# Si el modelo se ajustó correctamente, almacenar los resultadosif (!is.null(model)) {# Extraer el BIC del modelo bic <-BIC(model)# Almacenar la información en la lista de resultados results[[length(results) +1]] <-list(predictors =paste(terms, collapse =" + "),model = model,BIC = bic ) } } }# Convertir la lista de resultados en un dataframe results_df <- results %>% purrr::map_df(function(res) {data.frame(predictors = res$predictors,BIC = res$BIC,stringsAsFactors =FALSE ) })# Ordenar los modelos por BIC de menor a mayor results_df <- results_df %>%arrange(BIC)return(results_df)}best_subset_multinom_interactions_parallel <-function(y, x.vars, data) {# y Nombre de la variable dependiente (cadena de texto)# x.vars Vector de nombres de predictores (caracter)# data Dataframe con los datos de entrenamiento# Cargar las librerías necesarias dentro de la funciónrequire(dplyr)require(purrr)require(tidyr)require(nnet)require(MASS)require(foreach)require(doParallel)require(progressr)# Iniciar los gestores de progresohandlers(global =TRUE)handlers("txt")# Generar todas las combinaciones posibles de predictores (efectos principales) main_effects_list <-lapply(1:length(x.vars), function(i) {combn(x.vars, i, simplify =FALSE) }) %>%unlist(recursive =FALSE)# Inicializar una lista para almacenar las fórmulas de los modelos formulas_list <-list()# Generar todas las fórmulas posibles con interacciones hasta de 3 variablesfor (main_effects in main_effects_list) {# Generar términos de interacción de hasta 3 variables interaction_terms <-character(0) # Aseguramos que es un vector de caracteres# Para interacciones de 2 variablesif (length(main_effects) >=2) { interaction_terms_2way <-combn(main_effects, 2, function(x) paste(x, collapse =":"), simplify =TRUE) interaction_terms <-c(interaction_terms, interaction_terms_2way) }# Para interacciones de 3 variablesif (length(main_effects) >=3) { interaction_terms_3way <-combn(main_effects, 3, function(x) paste(x, collapse =":"), simplify =TRUE) interaction_terms <-c(interaction_terms, interaction_terms_3way) }# Generar todas las combinaciones posibles de efectos principales main_effects_subsets <-lapply(1:length(main_effects), function(i) {combn(main_effects, i, simplify =FALSE) }) %>%unlist(recursive =FALSE)# Para cada subconjunto de efectos principalesfor (me in main_effects_subsets) {# Iniciar con los efectos principales terms <- me# Identificar interacciones cuyos efectos principales están en 'me'if (length(interaction_terms) >0) { possible_interactions <- interaction_terms[vapply(interaction_terms, function(x) { vars_in_interaction <-unlist(strsplit(x, ":"))all(vars_in_interaction %in% me) }, FUN.VALUE =logical(1)) ] } else { possible_interactions <-character(0) }# Generar todas las combinaciones posibles de estas interacciones interaction_subsets <-list(character(0)) # Incluir el caso sin interaccionesif (length(possible_interactions) >0) { interaction_combinations <-lapply(1:length(possible_interactions), function(i) {combn(possible_interactions, i, simplify =FALSE) }) %>%unlist(recursive =FALSE) interaction_subsets <-c(interaction_subsets, interaction_combinations) }# Para cada combinación de interaccionesfor (ints in interaction_subsets) { full_terms <-c(terms, ints)# Crear la fórmula del modelo y almacenarla formula_str <-paste(y, "~", paste(full_terms, collapse ="+")) formulas_list <-append(formulas_list, list(formula_str)) } } }# Eliminar posibles duplicados de fórmulas formulas_list <-unique(formulas_list)# Total de modelos a ajustar total_models <-length(formulas_list)# Iniciar el progreso p <-progressor(steps = total_models)# Ajustar los modelos en paralelo usando foreach results_list <-foreach(i =1:total_models, .packages =c("nnet", "MASS"), .combine ='rbind') %dopar% { formula_str <- formulas_list[[i]] formula <-as.formula(formula_str)# Ajustar el modelo model <-tryCatch( nnet::multinom(formula, data = data, trace =FALSE),error =function(e) NULL,warning =function(w) NULL )# Actualizar el progresop(sprintf("Ajustando modelo %d de %d", i, total_models))# Si el modelo se ajustó correctamente, almacenar los resultadosif (!is.null(model)) { bic <-BIC(model)data.frame(predictors = formula_str,BIC = bic,stringsAsFactors =FALSE ) } else {NULL } }# Convertir los resultados a dataframe y ordenar por BIC results_df <-as.data.frame(results_list) results_df <- results_df %>%arrange(BIC)return(results_df)}num_cores <- parallel::detectCores() -1cl <-makeCluster(num_cores)registerDoParallel(cl)#pacman job kableExtra tidyverse cluster WeightedCluster devtools TraMineR TraMineRextras NbClust haven ggseqplot gridExtra Tmisc factoextra reticulate withr rmarkdown quartooptions(knitr.kable.NA ='')#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:##:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#knitr::knit_hooks$set(time_it =local({ now <-NULLfunction(before, options) {if (before) {# record the current time before each chunk now <<-Sys.time() } else {# calculate the time difference after a chunk res <-ifelse(difftime(Sys.time(), now)>(60^2),difftime(Sys.time(), now)/(60^2),difftime(Sys.time(), now)/(60^1))# return a character string to show the time x<-ifelse(difftime(Sys.time(), now)>(60^2),paste("Tiempo que demora esta sección:", round(res,1), "horas"),paste("Tiempo que demora esta sección:", round(res,1), "minutos"))paste('<div class="message">', gsub('##', '\n', x),'</div>', sep ='\n') } }}))knitr::opts_chunk$set(time_it =TRUE)#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:format_cells <-function(df, rows ,cols, value =c("italics", "bold", "strikethrough")){# select the correct markup# one * for italics, two ** for bold map <-setNames(c("*", "**", "~~"), c("italics", "bold", "strikethrough")) markup <- map[value] for (r in rows){for(c in cols){# Make sure values are not factors df[[c]] <-as.character( df[[c]])# Update formatting df[r, c] <-ifelse(nchar(df[r, c])==0,"",paste0(markup, gsub(" ", "", df[r, c]), markup)) } }return(df)}#To produce line breaks in messages and warningsknitr::knit_hooks$set(error =function(x, options) {paste('\n\n<div class="alert alert-danger">',gsub('##', '\n', gsub('^##\ Error', '**Error**', x)),'</div>', sep ='\n') },warning =function(x, options) {paste('\n\n<div class="alert alert-warning">',gsub('##', '\n', gsub('^##\ Warning:', '**Warning**', x)),'</div>', sep ='\n') },message =function(x, options) {paste('<div class="message">',gsub('##', '\n', x),'</div>', sep ='\n') })#_#_#_#_#_#_#_#_#_#_#_#_#_invisible("Function to format CreateTableOne into a database")as.data.frame.TableOne <-function(x, ...) {capture.output(print(x,showAllLevels =TRUE, varLabels = T,...) -> x) y <-as.data.frame(x) y$characteristic <- dplyr::na_if(rownames(x), "") y <- y %>%fill(characteristic, .direction ="down") %>% dplyr::select(characteristic, everything())rownames(y) <-NULL y}#_#_#_#_#_#_#_#_#_#_#_#_#_# Austin, P. C. (2009). The Relative Ability of Different Propensity # Score Methods to Balance Measured Covariates Between # Treated and Untreated Subjects in Observational Studies. Medical # Decision Making. https://doi.org/10.1177/0272989X09341755smd_bin <-function(x,y){ z <- x*(1-x) t <- y*(1-y) k <-sum(z,t) l <- k/2return((x-y)/sqrt(l))}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:if(.Platform$OS.type =="windows") withAutoprint({memory.size()memory.size(TRUE)memory.limit()})
invisible("info de validación de modelos con bootstrap")# PAM Trimestral= La solución de 4 y 6 o 7 cluster parecen tener buenos índices de calidad. # De todas formas, los índices ASW se encuentran en niveles que reflejan buena calidad, # 6 - 7 cluster reflejan de mejor forma las distancias entre los puntos.# # Comb= ASW todos sobre umbral (density, todos sobre distribución); HC todos sobre # umbral (density, 7 a 10, 3 y 2 sobre el umbral); HG debajo umbral (density, 3 a 15 # debajo distribución, 2 en distribución); PBC debajo umbral (density, 12 y 13 en distribución# , 14 y 15 sobre, 2 a 11, debajo)# Seq= ASW todos sobre umbral (2 en distribución, 3 en adelante, sobre distribución); # HC 2 a 9 sobre umbral (density, todos debajo); 6 en adelante sobre umbral HG # (density, 5 en adelante sobre distribución, resto en distribución); PBC sobre umbral # desde 5 en adelante (density, sobre distribución desde 6 en adelante) # 6 EN ADELANTE PODRÍA SER DISTINTO A NULO# invisible("Hacemos clasificación de pertenencia cluster y ponemos etiquetas")ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam <-factor(pamRange_quarter_om$clustering$cluster6,levels=rev(attr( sort(table(pamRange_quarter_om$clustering$cluster6)), "name")),labels=c("6035_Un evento, TSM", "6025_Un evento, TUS", "5939_Un evento TSM larga duración", "5989_Un evento, comorbilidad", "6036_TSM, 1 año después, otras causas", "5710_TSM, 1 año después, TSM"))
Tiempo que demora esta sección: 0 minutos
Vemos los diagnósticos que vienen después de aquellos cluster con más de un ingreso.
Code
diag_6036<-df_filled %>% dplyr::filter(run %in%subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, clus_pam=="6036_TSM, 1 año después, otras causas")$run) %>% dplyr::select(run, diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, diag10, diag11, fecha_egreso_rec_fmt, estab_homo) %>% dplyr::group_by(run) %>% dplyr::filter(row_number() !=1) %>%# Elimina la primera observación de cada run dplyr::mutate(all_diags =paste(na.omit(c(diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, diag10, diag11)), collapse =", ") ) %>% dplyr::summarise(all_diags =first(all_diags),fecha_egreso_rec_fmt =first(fecha_egreso_rec_fmt),estab_homo =first(estab_homo) ) %>% dplyr::ungroup() %>% dplyr::pull(all_diags) %>%# Extraer la columna all_diags como vectorstrsplit(split =", ") %>%# Separar cada diagnóstico por comasunlist() # Convertir la lista en un vector# 1. **F192** (25, 2.96%): Trastorno mental y del comportamiento debido al uso de múltiples drogas y al uso de otras sustancias psicoactivas, síndrome de abstinencia.# 2. **F609** (17, 2.01%): Trastorno específico de la personalidad no especificado.# 3. **O800** (17, 2.01%): Parto espontáneo por vía vaginal.# 4. **F603** (16, 1.89%): Trastorno de la personalidad emocionalmente inestable, tipo impulsivo.# 5. **E101** (15, 1.77%): Diabetes mellitus insulino-dependiente, no controlada.# 6. **Z370** (14, 1.65%): Nacimiento de un solo feto nacido vivo.# 7. **F322** (13, 1.54%): Episodio depresivo grave sin síntomas psicóticos.# 8. **G409** (13, 1.54%): Epilepsia, no especificada.# 9. **N390** (13, 1.54%): Infección del tracto urinario, sitio no especificado.# 10. **Z518** (13, 1.54%): Otras atenciones médicas especificadas.# 11. **T509** (10, 1.18%): Envenenamiento por otros psicotrópicos y sustancias psicotrópicas no especificadas.# 12. **G629** (8, 0.95%): Polineuropatía, no especificada.# 13. **K358** (8, 0.95%): Apendicitis crónica.# 14. **N10X** (8, 0.95%): Nefritis tubulointersticial aguda.# 15. **O829** (8, 0.95%): Parto por cesárea no especificado.# 16. **F191** (7, 0.83%): Trastorno mental y del comportamiento debido al uso de múltiples drogas y al uso de otras sustancias psicoactivas, dependencia actual.# 17. **F329** (7, 0.83%): Episodio depresivo no especificado.# 18. **Z291** (7, 0.83%): Asesoramiento sobre la prevención del abuso de sustancias psicoactivas.# 19. **F199** (6, 0.71%): Trastorno mental y del comportamiento debido al uso de múltiples drogas y al uso de otras sustancias psicoactivas, no especificado.# 20. **K808** (6, 0.71%): Cálculo de la vesícula biliar sin colecistitis.# 21. **O470** (6, 0.71%): Falso trabajo de parto.# 22. **O809** (6, 0.71%): Parto por vía vaginal no especificado.# 23. **O821** (5, 0.59%): Parto por cesárea electiva.# 24. **R458** (5, 0.59%): Otros síntomas y signos que involucran el estado emocional.# 25. **X590** (5, 0.59%): Exposición a otros factores no especificados.# 26. **Z910** (5, 0.59%): Antecedentes personales de riesgo no especificado.# invisible("6035 (n=4476): trayectoria con sólo un evento hosp por TSM")invisible("6025 (n=680): trayectoria con sólo un evento hosp. con TUS")invisible("5939 (n=319): trayectoria con 2 trimestres continuis por morbilidad psiquiátrica")invisible("5989 (n=207): trayectoria de sólo un evento en por comorbilidad")invisible("6036 (n=202): trayectoria de morbilidad, pero a los 4 trimestres posteriores tienen ingresos por otras causas")invisible("cuáles son esas causas principalmente??, son distintas de las que tiene la población?")invisible("5710 (n=154): trayectoria de morbilidad, pero a los 4 trimestres posteriores tienen ingresos por morbilidad")invisible("qué morbilidad principalmente??, son distintas de las que tiene la población?")
Tiempo que demora esta sección: 0 minutos
Code
diag_5710<- df_filled %>% dplyr::filter(run %in%subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, clus_pam=="5710_TSM, 1 año después, TSM")$run) %>% dplyr::select(run, diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, diag10, diag11, fecha_egreso_rec_fmt, estab_homo) %>% dplyr::group_by(run) %>% dplyr::filter(row_number() !=1) %>%# Elimina la primera observación de cada run dplyr::mutate(all_diags =paste(na.omit(c(diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, diag10, diag11)), collapse =", ") ) %>% dplyr::summarise(all_diags =first(all_diags),fecha_egreso_rec_fmt =first(fecha_egreso_rec_fmt),estab_homo =first(estab_homo) ) %>% dplyr::ungroup() %>% dplyr::pull(all_diags) %>%# Extraer la columna all_diags como vectorstrsplit(split =", ") %>%# Separar cada diagnóstico por comasunlist() # Convertir la lista en un vector# # 1. **F603** (71, 5.97%): Trastorno de la personalidad emocionalmente inestable, tipo impulsivo.# 2. **F609** (54, 4.54%): Trastorno específico de la personalidad no especificado.# 3. **F329** (51, 4.29%): Episodio depresivo no especificado.# 4. **F322** (47, 3.95%): Episodio depresivo grave sin síntomas psicóticos.# 5. **F319** (43, 3.61%): Trastorno depresivo recurrente, episodio no especificado.# 6. **F209** (42, 3.53%): Esquizofrenia, no especificada.# 7. **F200** (40, 3.36%): Esquizofrenia paranoide.# 8. **F192** (33, 2.77%): Trastorno mental y del comportamiento debido al uso de múltiples drogas y al uso de otras sustancias psicoactivas, síndrome de abstinencia.# 9. **C490** (29, 2.44%): Tumor maligno de tejido mesotelial y tejido blando, sitio no especificado.# 10. **G909** (21, 1.76%): Trastorno no especificado del sistema nervioso autónomo.# 11. **Z511** (20, 1.68%): Atención para tratamiento de tumores malignos.# 12. **F432** (17, 1.43%): Reacción de estrés agudo.# 13. **F431** (16, 1.34%): Trastorno de estrés postraumático.# 14. **F449** (15, 1.26%): Trastorno de ansiedad no especificado.# 15. **F070** (14, 1.18%): Trastornos orgánicos de la personalidad.# 16. **Z915** (14, 1.18%): Historia personal de lesión autoinfligida.# 17. **F191** (12, 1.01%): Trastorno mental y del comportamiento debido al uso de múltiples drogas y al uso de otras sustancias psicoactivas, dependencia actual.# 18. **E669** (11, 0.92%): Obesidad no especificada.# 19. **T424** (11, 0.92%): Intoxicación por otros psicotrópicos.# 20. **F608** (10, 0.84%): Otros trastornos específicos de la personalidad.# 21. **G409** (10, 0.84%): Epilepsia, no especificada.
Generamos un gráfico de PPOO por cada conglomerado.
Code
ppoo_clus<- df_filled[,c("run","glosa_pueblo_originario")] %>% dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run", "clus_pam","factor_inclusivo_real_hist_mas_autperc")], by="run", multiple="first") %>% dplyr::mutate(glosa_pueblo_originario_rec= dplyr::case_when(glosa_pueblo_originario=="NINGUNO"& factor_inclusivo_real_hist_mas_autperc!="00"~"DESCONOCIDO", T~glosa_pueblo_originario)) %>% janitor::tabyl(glosa_pueblo_originario_rec, clus_pam) %>% janitor::adorn_percentages("row")reshape2::melt(ppoo_clus, id.vars ="glosa_pueblo_originario_rec") %>% dplyr::mutate(glosa_pueblo_originario_rec= dplyr::recode(glosa_pueblo_originario_rec, "OTRO (ESPECIFICAR)"="OTRO(n=77)", "RAPA NUI (PASCUENSE)"="RAPA NUI(n=34)", "YAGÁN (YÁMANA)"="YAGÁN(n=2)","AYMARA"="AYMARA(n=13)","COLLA"="COLLA(n=6)","DIAGUITA"="DIAGUITA(n=3)","KAWÉSQAR"="KAWÉSQAR(n=4)","MAPUCHE"="MAPUCHE(n=255)","DESCONOCIDO"=".DESCONOCIDO(n=1.985)","NINGUNO"=".NINGUNO(n=9.156)")) %>%#dplyr::filter(glosa_pueblo_originario!="NINGUNO") %>% # dplyr::mutate(variable= # dplyr::recode(variable, "6036_TSM, 1 año después, otras causas"="6036_TSM, 1 año\ndespués, otras causas",# "5710_TSM, 1 año después, TSM"="5710_TSM, 1 año\ndespués, TSM",# "5939_Un evento TSM larga duración"="5939_Un evento TSM\nlarga duración",# "5989_Un evento, comorbilidad"="5989_Un evento,\ncomorbilidad")) %>% ggplot(aes(x = glosa_pueblo_originario_rec, y = value, fill = variable)) +geom_bar(stat ="identity", position ="fill") +scale_fill_manual(values =c("#E27A5B","#6B8E23", "#D2B48C", "#696969", "#BDB76B", "#4682B4")) +labs(title =NULL,x ="Grupo Étnico",y ="Proporción de Casos",fill ="Grupos") +# Cambia el título de la leyenda a "Grupos"theme_minimal() +theme(axis.text.y =element_text(size =12), # Tamaño de las etiquetas de los grupos étnicosaxis.text.x =element_text(size =12), # Tamaño de las etiquetas del eje Xaxis.title.x =element_text(size =14), # Tamaño del título del eje Xaxis.title.y =element_text(size =14), # Tamaño del título del eje Yplot.title =NULL, # Tamaño y estilo del título del gráficolegend.title =element_text(size =14, margin =margin(b =-.1)), # Tamaño del título de la leyendalegend.spacing.y =unit(1.5, "lines"),legend.box.spacing =unit(0.5, "lines"), # Controla el espacio entre la leyenda y el gráficolegend.margin =margin(5, 5, 5, 5), legend.key.height =unit(1, "cm"), legend.text =element_text(size =12) # Tamaño del texto de la leyenda ) +coord_flip() # Hacer el gráfico horizontalggsave("grafico_ancho_achatado.png", width =8, height =3, dpi=1000)
PPOO por cluster
Tiempo que demora esta sección: 0 minutos
Vemos los gráficos de las trayectorias
Code
categories<-attr(States_Wide.seq_quarter_t_prim_adm_cens, "labels")new_labels <- categoriesnew_labels[which(categories =="Otras causas")] <-"Otras\ncausas"#new_labels[which(categories == "Consumo\nde sustancias")] <- "Consumo de\nsustancias"seq_plot <-ggseqiplot(States_Wide.seq_quarter_t_prim_adm_cens, group=factor(pamRange_quarter_om$clustering$cluster6,levels=rev(attr( sort(table(pamRange_quarter_om$clustering$cluster6)), "name"))), facet_ncol=3, facet_nrow=3) +theme(legend.position ="none")+labs(x="Trimestres", y="# IDs de usuarios")+#guides(fill = guide_legend(nrow = 1))+theme(axis.text.y =element_text(size =15), # Tamaño de las etiquetas de los grupos étnicosaxis.text.x =element_text(size =15), # Tamaño de las etiquetas del eje Xaxis.title.x =element_text(size =15), # Tamaño del título del eje Xaxis.title.y =element_text(size =15, margin =margin(r =-10)),#,margin = margin(l = -10)),strip.text =element_text(size =16, margin =margin(b =-15)),legend.text =element_text(size =15),legend.spacing.x =unit(0.1, 'cm'), # Alinea el título de la leyenda hacia la izquierdalegend.box.margin =margin(t =0, r =0, b =0, l =-50),legend.position ="bottom", legend.justification ="left",panel.spacing.y =unit(0.5, "lines")#legend.key.size = unit(1.5, "lines"), # Aumenta el tamaño de los símbolos en la leyenda )+guides(fill =guide_legend(nrow =1)) +scale_fill_manual(labels = new_labels, values=c("#E2725B", "#556B2F", "#D2B48C",#"#8B4513","#FFFFFF","#808080","#000000"))seq_plotggsave(filename="clusterspam_om6_mod.png",seq_plot, width =8.5, height =5.5, dpi=1000)
Trayectorias de hospitalización, orden de sujetos según el primer estado observado y su duración, representando a cada individuo como una línea en el gráfico.
Tiempo que demora esta sección: 0.4 minutos
Code
seq_plot2 <-ggseqdplot(States_Wide.seq_quarter_t_prim_adm_cens, group=factor(pamRange_quarter_om$clustering$cluster6,levels=rev(attr( sort(table(pamRange_quarter_om$clustering$cluster6)), "name"))), facet_ncol=3, facet_nrow=3) +theme(legend.position ="none")+# Colocar la leyenda abajolabs(x="Trimestres", y="Frecuencia relativa de estados")+theme(axis.text.y =element_text(size =15), # Tamaño de las etiquetas de los grupos étnicosaxis.text.x =element_text(size =15), # Tamaño de las etiquetas del eje Xaxis.title.x =element_text(size =15), # Tamaño del título del eje Xaxis.title.y =element_text(size =15, margin =margin(r =-5)),strip.text =element_text(size =15),panel.spacing.y =unit(0.5, "lines") ) # Colocar la leyenda abajoseq_plot2ggsave("clusterspam_om62_mod.png",seq_plot2, width =8, height =5.5, dpi=1000)
Trayectorias de hospitalización, frecuencia relativa de estados en un gráfico de barras apiladas por trimestre.
Tiempo que demora esta sección: 0.1 minutos
De este modo, presenta el cambio agregado en la distribución de estados a lo largo del tiempo, sin considerar las secuencias individuales.
2.1.1.Exploración transiciones
2.1.1.a Transiciones- RM y no RM
Tasas de transición no RM a RM y viceversa
Code
invisible("Tasas de transición no RM a RM y viceversa")trim_tasa_pam_clus6_cens_cnt<-seqcount_t(States_Wide.seq_quarter_t_prim_adm_RM_cens, group=factor(pamRange_quarter_om$clustering$cluster6,levels=rev(attr( sort(table(pamRange_quarter_om$clustering$cluster6)), "name")))) %>% dplyr::filter(count>0) %>% dplyr::mutate(trans =paste0(from,"_", to)) %>% dplyr::mutate(across(c("from","to"),~gsub("\\[->\\s*|\\s*->\\s*\\]|\\[|\\]", "", .))) trim_tasa_pam_clus6_cens_rate<-seqtrate_t(States_Wide.seq_quarter_t_prim_adm_RM_cens, group=factor(pamRange_quarter_om$clustering$cluster6,levels=rev(attr( sort(table(pamRange_quarter_om$clustering$cluster6)), "name")))) %>% dplyr::filter(rate>0) %>% dplyr::mutate(trans =paste0(from,"_", to)) %>% dplyr::mutate(across(c("from","to"),~gsub("\\[->\\s*|\\s*->\\s*\\]|\\[|\\]", "", .)))
Tiempo que demora esta sección: 0 minutos
Code
trim_tasa_pam_clus6_cens_rate %>% dplyr::left_join(trim_tasa_pam_clus6_cens_cnt, by=c("from"="from", "glosa_sexo"="glosa_sexo","to"="to")) %>% dplyr::rename("recuento"="count") %>% dplyr::filter(from %in%c("RM", "noRM")) %>% dplyr::mutate(glosa_sexo=factor(glosa_sexo,levels=rev(attr( sort(table(pamRange_quarter_om$clustering$cluster6)), "name")))) %>%ggplot(aes(x = from, y = to, fill = rate, size=log(recuento+1))) +geom_tile() +coord_flip()+scale_fill_gradient(low ="white", high ="blue") +# Ajusta la escala de colores según tus preferenciaslabs(title ="Tasas de transición, Trimestre (s/censura)",x ="Desde",y ="Hacia",fill ="Rate") +theme_minimal() +facet_wrap(~glosa_sexo)+theme(axis.text.x =element_text(angle =45, hjust =1))+geom_text(aes(label =sprintf("%1.2f", rate), size =log(recuento+1)*.5), color ="black")invisible("Hay muy pocos casos que se entrecruzan entre noRM y RM (fuera de la diagnonal)")
Porcentajes de transición no-RM y RM por cada cluster
Tiempo que demora esta sección: 0.1 minutos
Hay muy pocos casos que se entrecruzan entre noRM y RM (fuera de la diagnonal)
trim_tasa2_pam_clus6_cens_rate %>% dplyr::left_join(trim_tasa2_pam_clus6_cens_cnt, by=c("from"="from", "glosa_sexo"="glosa_sexo","to"="to")) %>% dplyr::rename("recuento"="count") %>%#dplyr::filter(from %in% c("RM", "noRM")) %>% dplyr::mutate(glosa_sexo=factor(glosa_sexo,levels=rev(attr( sort(table(pamRange_quarter_om$clustering$cluster6)), "name")))) %>%ggplot(aes(x = from, y = to, fill = rate, size=log(recuento+1))) +geom_tile() +coord_flip()+scale_fill_gradient(low ="white", high ="blue") +# Ajusta la escala de colores según tus preferenciaslabs(title ="Tasas de transición, Trimestre (s/censura)",x ="Desde",y ="Hacia",fill ="Rate") +theme_minimal() +facet_wrap(~glosa_sexo)+theme(axis.text.x =element_text(angle =45, hjust =1))+geom_text(aes(label =sprintf("%1.2f", rate), size =log(recuento+1)*.5), color ="black")
Porcentajes de transición no-RM y RM por cada cluster
Tiempo que demora esta sección: 0.3 minutos
2.1.1.c Tiempo promedio por cluster
Code
seq_mean_t(States_Wide.seq_quarter_t_prim_adm_cens, factor(pamRange_quarter_om$clustering$cluster6,levels=rev(attr( sort(table(pamRange_quarter_om$clustering$cluster6)), "name"))))%>% data.table::as.data.table(keep.rowname=T) %>%ggplot(aes(x=rn, fill= factor_inclusivo, y=Mean))+geom_bar(width =1, stat ="identity") +theme_minimal() +labs(title =NULL,x =NULL,y =NULL) +coord_flip()+theme(#axis.text.x = element_blank(),#axis.text.y = element_blank(),panel.grid =element_blank()) +# scale_fill_brewer(palette = "Pastel1", labels=c("Sin\nautoidentificación\nni reconocimiento", "Autoidentificación\nsin reconocimiento", "Ambas")) +geom_text(aes(label =round(Mean,1)), position =position_stack(vjust =0.5), size =2.5, # Ajusta el tamaño de la fuente aquícolor ="black", # Color del textofamily ="sans", # Puedes cambiar la fuente si lo deseasbackground =element_rect(fill ="white", color =NA)) +# Fondo blancotheme(legend.title =element_blank())invisible("No me aporta mucho")# seq_mean_t_dos_grupos(States_Wide.seq_quarter_t_prim_adm_cens, group1=ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam, group2=ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$factor_inclusivo_real_hist_mas_autperc)
Tiempo promedio en cada estado por estatus PPOO (Trimestral c/censura)
Tiempo que demora esta sección: 0.1 minutos
2.1.2. Propiedades secuencias
Vemos el número de transiciones, número de subsecuencias y la entropía
Code
## number of transitions ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$ntrans <-as.numeric(seqtransn(States_Wide.seq_quarter_t_prim_adm))## number of subsequencesing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$subsn <-as.numeric(seqsubsn(States_Wide.seq_quarter_t_prim_adm))ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$entropy_long <-as.numeric(seqient(States_Wide.seq_quarter_t_prim_adm))ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens %>% dplyr::group_by(clus_pam) %>% dplyr::summarise(promedio_entropia=mean(entropy_long), mediana_entropia=quantile(entropy_long, .5),p25_entropia=quantile(entropy_long, .25),p75_entropia=quantile(entropy_long, .75),promedio_n_trans=mean(ntrans), #Las transiciones se refieren a los cambios de estado en las secuencias.sd_n_trans=sd(ntrans), #subsecuencia es una parte de una secuencia más larga que respeta el orden temporal original, #pero no necesariamente contiene todos los estadospromedio_subsn=mean(subsn), sd_subsn=sd(subsn)) %>% dplyr::mutate_if(is.numeric, ~sprintf("%1.2f", .)) %>% knitr::kable("html", col.names =c("Cluster", "Promedio Entropía", "Mediana Entropía", "Entropía Percentil 25", "Entropía Percentil 75", "Promedio Número de Transiciones", "Desviación Estándar Número de Transiciones", "Promedio Número de Subsecuencias", "Desviación Estándar Número de Subsecuencias"), caption="Resumen de Métricas por Cluster")
Resumen de Métricas por Cluster
Cluster
Promedio Entropía
Mediana Entropía
Entropía Percentil 25
Entropía Percentil 75
Promedio Número de Transiciones
Desviación Estándar Número de Transiciones
Promedio Número de Subsecuencias
Desviación Estándar Número de Subsecuencias
6035_Un evento, TSM
0.16
0.12
0.12
0.20
1.71
1.20
9.38
46.91
6025_Un evento, TUS
0.18
0.12
0.12
0.25
1.78
1.34
10.95
27.85
5939_Un evento TSM larga duración
0.28
0.26
0.20
0.32
2.24
1.68
17.64
52.58
5989_Un evento, comorbilidad
0.20
0.12
0.12
0.25
2.00
1.39
13.92
41.49
6036_TSM, 1 año después, otras causas
0.33
0.32
0.25
0.38
4.35
1.63
50.74
76.27
5710_TSM, 1 año después, TSM
0.37
0.38
0.27
0.43
4.62
1.99
57.57
83.32
Tiempo que demora esta sección: 0 minutos
2.1.3. Comparación variables
2.1.3.a. Comparación covariables- PPOO
Code
round(prop.table(table(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam, ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$factor_inclusivo_real_hist_mas_autperc),1),2)ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens %>%count(clus_pam, factor_inclusivo_real_hist_mas_autperc) %>%group_by(clus_pam) %>%mutate(prop = scales::percent(n /sum(n))) %>%select(-n) %>%pivot_wider(names_from = factor_inclusivo_real_hist_mas_autperc, values_from = prop, values_fill ="0") %>% knitr::kable("markdown", col.names=c("Conglomerados","No se identifica/no pertenece", "No se identifica/hay reconocimiento", "Se identifica/hay reconocimeinto"), caption="Porcentajes por fila, conglomerado vs. Pertenencia/identificación + Reconocimento CONADI PPOO")# 00 10 11# 6035 0.81 0.11 0.08# 6025 0.79 0.12 0.10# 5939 0.83 0.11 0.07# 5989 0.84 0.09 0.08# 6036 0.80 0.11 0.09# 5710 0.79 0.12 0.09invisible("6025 tiwnw un poxo mas PPOO, lo mismo con 5710")
Pearson's Chi-squared test
data: .
X-squared = 5.6334, df = 10, p-value = 0.8451
Fisher's Exact Test for Count Data with simulated p-value (based on
1e+05 replicates)
data: .
p-value = 0.8465
alternative hypothesis: two.sided
round(prop.table(table(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam, ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$factor_inclusivo_real_hist_mas_autperc_bin),1),2) %>% knitr::kable("markdown", col.names=c("Conglomerados","No se identifica/no pertenece", "Se identifica/pertenece"), caption="Porcentajes por fila, conglomerado vs. Pertenencia/identificación PPOO")
Porcentajes por fila, conglomerado vs. Pertenencia/identificación PPOO
Conglomerados
No se identifica/no pertenece
Se identifica/pertenece
6035_Un evento, TSM
0.81
0.19
6025_Un evento, TUS
0.79
0.21
5939_Un evento TSM larga duración
0.83
0.17
5989_Un evento, comorbilidad
0.84
0.16
6036_TSM, 1 año después, otras causas
0.80
0.20
5710_TSM, 1 año después, TSM
0.79
0.21
Tiempo que demora esta sección: 0 minutos
Hicimos una prueba post-hoc usando Bonferroni
Code
# # Beasley, T. M., & Schumacker, R. E. (1995). Multiple Regression Approach to Analyzing Contingency# Tables: Post Hoc and Planned Comparison Procedures. The Journal of Experimental Education, 64(1), 79–93. https://doi.org/10.1080/00220973.1995.9943797chisq.posthoc.test(table(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam, ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$factor_inclusivo_real_hist_mas_autperc_bin), method ="bonferroni") %>%pivot_longer(cols =c(`0`, `1`), names_to ="Group", values_to ="Values") %>%pivot_wider(names_from = Value, values_from = Values) %>%separate(Dimension, into =c("Code", "Evento"), sep ="_", extra ="merge") %>% dplyr::filter(Group==1) %>% dplyr::mutate(Residuals=sprintf("%1.2f", Residuals)) %>%select(Code, Evento, Residuals, `p values`) %>% knitr::kable("markdown", col.names=c("Código","Descripción", "Residuos", "Sig."), caption="Post-hoc, conglomerado vs. Pertenencia/identificación PPOO")
Post-hoc, conglomerado vs. Pertenencia/identificación PPOO
Código
Descripción
Residuos
Sig.
6035
Un evento, TSM
-0.30
1
6025
Un evento, TUS
1.41
1
5939
Un evento TSM larga duración
-1.03
1
5989
Un evento, comorbilidad
-1.12
1
6036
TSM, 1 año después, otras causas
0.31
1
5710
TSM, 1 año después, TSM
0.42
1
Tiempo que demora esta sección: 0 minutos
2.1.3.b. Comparación covariables- Mortalidad
Code
# invisible("No hay nada, el tiempo promedio de censura es similar")ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens %>% dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)) %>% janitor::tabyl(clus_pam,death_time_rec) %>% janitor::adorn_percentages("row")%>% dplyr::mutate(`1`=scales::percent(`1`, accuracy=.1)) %>% dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens %>% dplyr::group_by(clus_pam) %>% dplyr::summarise(mean=sprintf("%1.1f",mean(cens_time))), by="clus_pam") %>% dplyr::select(-`0`) %>% knitr::kable("markdown", col.names=c("Conglomerado","Mortalidad observada", "Promedio"), caption="Post-hoc, conglomerado vs. Pertenencia/identificación PPOO")
Post-hoc, conglomerado vs. Pertenencia/identificación PPOO
Conglomerado
Mortalidad observada
Promedio
6035_Un evento, TSM
0.9%
17.9
6025_Un evento, TUS
1.9%
18.2
5939_Un evento TSM larga duración
2.2%
18.0
5989_Un evento, comorbilidad
1.9%
18.1
6036_TSM, 1 año después, otras causas
1.0%
17.8
5710_TSM, 1 año después, TSM
1.9%
18.1
Tiempo que demora esta sección: 0 minutos
Code
ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens %>% dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)) %>% janitor::tabyl(death_time_rec,clus_pam) %>% janitor::chisq.test(correct=T)#X-squared = 10.621, df = 5, p-value = 0.05943ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens %>% dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)) %>% janitor::tabyl(death_time_rec,clus_pam) %>% janitor::fisher.test(simulate.p.value=T, B=1e5)#p-value = 0.03169invisible("no se basa en la distribución chi-cuadrado. Fisher se basa en permutaciones exactas, por lo que no se calculan df.")invisible("Podría haber algo aquí, aunque son números pequeños")invisible("6036 (morbilidad y a los 4 meses otras causas) y 6035 (sólo un evento por trno.SM) tienen un 1% de gente que muere antes")invisible("5939 (dos semestres continuos por morbilidad psiquiátrica) en cambio, tiene 2.2%")#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_##_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_tabla_cluster_mortalidad<- ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens %>% dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1))chisq.posthoc.test(table(tabla_cluster_mortalidad$clus_pam, tabla_cluster_mortalidad$death_time_rec))%>%pivot_longer(cols =c(`0`, `1`), names_to ="Group", values_to ="Values") %>%pivot_wider(names_from = Value, values_from = Values) %>%separate(Dimension, into =c("Code", "Evento"), sep ="_", extra ="merge") %>% dplyr::filter(Group==1) %>% dplyr::mutate(Residuals=sprintf("%1.2f", as.numeric(Residuals))) %>%select(Code, Evento, Residuals, `p values`) %>% knitr::kable("markdown", col.names=c("Código","Descripción", "Residuos", "Sig."), caption="Post-hoc, conglomerado vs. mortalidad")
Pearson's Chi-squared test
data: .
X-squared = 10.621, df = 5, p-value = 0.05943
Fisher's Exact Test for Count Data with simulated p-value (based on
1e+05 replicates)
data: .
p-value = 0.03191
alternative hypothesis: two.sided
Post-hoc, conglomerado vs. mortalidad
Código
Descripción
Residuos
Sig.
6035
Un evento, TSM
-2.99
0.0335*
6025
Un evento, TUS
1.95
0.6201
5939
Un evento TSM larga duración
1.77
0.9118
5989
Un evento, comorbilidad
1.06
1
6036
TSM, 1 año después, otras causas
-0.23
1
5710
TSM, 1 año después, TSM
0.93
1
Tiempo que demora esta sección: 0 minutos
2.1.3.c. Comparación covariables- no RM vs. RM
Code
round(prop.table(table(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam, ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$codigo_region_rec_base),1),2) %>% knitr::kable("markdown", caption="Porcentajes por fila")
Edad promedio primer ingreso con intervalo de confianza por conglomerado
Tiempo que demora esta sección: 0.1 minutos
Code
anova <-oneway.test(min_edad_anos ~ clus_pam, data = dt_ing_calendar_quarter_t_desde_primera_adm_dedup %>% dplyr::filter(quarter ==0) %>% dplyr::inner_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run","clus_pam")], by="run"),var.equal = F)# Ver los resultados del ANOVAprint(anova)
One-way analysis of means (not assuming equal variances)
data: min_edad_anos and clus_pam
F = 31.673, num df = 5.00, denom df = 604.01, p-value < 2.2e-16
Pearson's Chi-squared test
data: .
X-squared = 34.049, df = 20, p-value = 0.0258
Fisher's Exact Test for Count Data with simulated p-value (based on
1e+05 replicates)
data: .
p-value = 0.02305
alternative hypothesis: two.sided
Porcentajes por columna, conglomerado vs. Beneficios
Pearson's Chi-squared test
data: .
X-squared = 3.0869, df = 5, p-value = 0.6866
Fisher's Exact Test for Count Data with simulated p-value (based on
1e+05 replicates)
data: .
p-value = 0.6785
alternative hypothesis: two.sided
Pearson's Chi-squared test
data: .
X-squared = 4.0391, df = 5, p-value = 0.5438
Fisher's Exact Test for Count Data with simulated p-value (based on
1e+05 replicates)
data: .
p-value = 0.5523
alternative hypothesis: two.sided
#Definimos la base de datos que agrupa por observación y nos permite unirla con nuestros clusterdias_ttos_base<-data_long_establecimiento_2024_std %>% dplyr::group_by(run) %>% dplyr::summarise(n_ttos=n(), promedio_dias=mean(days_elapsed))invisible("Prueba de Levene par igualdad de varianzas")with(dplyr::inner_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, dias_ttos_base, by =c("run"="run"), multiple ="first"), car::leveneTest(n_ttos, clus_pam))# Realizar el ANOVA comparando la edad media entre los diferentes conglomerados (clus_pam)anova_n_ttos <-oneway.test(n_ttos ~ clus_pam, data = ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens%>% dplyr::inner_join(dias_ttos_base, by =c("run"="run"), multiple ="first"), var.equal = F)# Ver los resultados del ANOVAprint(anova_n_ttos)#F = 101.98, num df = 5.00, denom df = 570.85, p-value < 2.2e-16rstatix::games_howell_test(n_ttos ~ clus_pam, data =ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens%>% dplyr::inner_join(dias_ttos_base, by =c("run"="run"), multiple ="first")) %>% dplyr::select(-1) %>% dplyr::mutate(summary =sprintf("%.2f [%.2f, %.2f], p= %s",as.numeric(estimate), as.numeric(conf.low), as.numeric(conf.high), ifelse(p.adj <0.001, "<0.001", sprintf("%.3f", p.adj)))) %>% dplyr::select(!any_of(c("estimate","conf.low", "conf.high", "p.adj", "p.adj.signif"))) %>% knitr::kable("markdown", col.names=c("Conglomerado1","Conglomerado2", "Estimación"), caption="Post-hoc, conglomerado vs. N°s días de tratamiento")
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 5 95.597 < 2.2e-16 ***
6032
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
One-way analysis of means (not assuming equal variances)
data: n_ttos and clus_pam
F = 101.98, num df = 5.00, denom df = 570.85, p-value < 2.2e-16
Post-hoc, conglomerado vs. N°s días de tratamiento
# Definir los datos correctamentedata <-data.frame(Grupo =c('1', '2', '3', '4', '5', '6'),Macrozona_Norte =c('+', '-', NA, NA, NA, NA),Macrozona_Sur =c('-', '+', NA, NA, NA, NA),RM =c(NA, '-', NA, '+', NA, NA),Sexo_Mujeres =c('+', '-', NA, '-', '+', NA),Edad =c(NA, '+', '-', NA, NA, NA),N_ttos =c('-', '-', NA, NA, NA, '+'),Dias_en_tto =c(NA, '-', NA, NA, NA, '+'))# Asegurar que los nombres de las columnas sean válidos y no haya espacios en blancocolnames(data) <-c('Grupo', 'Macrozona\nNorte', 'Macrozona\nSur', 'RM', 'Sexo\n(Mujeres)', 'Edad', 'N°ttos', 'Dias_en\ntto.')# Derretir el dataframe para que sea adecuado para ggplot2data_melt <- reshape2::melt(data, id.vars ='Grupo', variable.name ='Variable', value.name ='Asociación')# Reemplazar los NA por un valor vacíodata_melt$Asociación[is.na(data_melt$Asociación)] <-"\n"# Crear el gráfico con ggplotdata_melt %>% dplyr::mutate(Variable =gsub("_", " ", Variable)) %>%ggplot(aes(x = Variable, y = Grupo, fill = Asociación)) +geom_tile(color ="white", size =0.8) +scale_fill_manual(values =c("+"="#556B2F", "-"="#E2725B", "\n"="white")) +labs(title =NULL, x ="Variables", y ="Conglomerado") +theme_minimal() +theme(#axis.text.x = element_text(angle = 45, hjust = 1),panel.grid =element_blank())+theme(axis.text.y =element_text(size =17, face ="bold"),#,margin = margin(l = 7)), # Tamaño de las etiquetas de los grupos étnicosaxis.text.x =element_text(size =17, face ="bold"), # Tamaño de las etiquetas del eje Xaxis.title.x =element_text(size =16, face ="bold"),#,margin = margin(t = -15)), # Tamaño del título del eje Xaxis.title.y =element_text(size =16, face ="bold"), # Tamaño del título del eje Yplot.title =NULL, # Tamaño y estilo del título del gráficolegend.title =element_text(size =17, face ="bold"), # Tamaño del título de la leyendalegend.spacing.y =unit(1.5, "lines"),legend.box.spacing =unit(0.5, "lines"), # Controla el espacio entre la leyenda y el gráficolegend.margin =margin(5, 5, 5, 5), legend.key.height =unit(1, "cm"), legend.text =element_text(size =15, face ="bold") # Tamaño del texto de la leyenda ) +coord_flip()ggsave("_figs/asociaciones.png", width=8.8*.8, height=5*.8, dpi=1000)
Número de tratamientos (promedio y bigotes en IC95%)
Tiempo que demora esta sección: 0.1 minutos
2.1.5. Regresión
Code
multinomial <- nnet::multinom( clus_pam ~ codigo_region_rec_base + glosa_sexo * factor_inclusivo_real_hist_mas_autperc_bin + prev_benef_rec_post,data = ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens %>% dplyr::inner_join(data_long_establecimiento_2024_std[,c("ESTAB_HOMO", "codigo_region", "nivel_de_atencion", "nivel_de_complejidad")], by =c("estab_homo_base"="ESTAB_HOMO"), multiple ="first") %>% dplyr::mutate(prev_benef_rec_post =fct_relevel(prev_benef_rec_post, "FONASA A")),maxit =1000, # Aumenta el número de iteracionestrace =TRUE)model_table <- gtsummary::tbl_regression(multinomial, exponentiate =TRUE)invisible("Usé la pertenencia a PPOO como binaria")model_table
# weights: 60 (45 variable)
initial value 10818.643675
iter 10 value 5847.307596
iter 20 value 5574.080348
iter 30 value 5525.522191
iter 40 value 5524.336545
iter 50 value 5524.304683
final value 5524.304601
converged
